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The fluctuation exchange, or FLEX, approximation for interacting electrons is applied to study 
instabilities in the standard three-band model for CUO2 layers in the high-temperature supercon- 
ductors. Both intra-orbital and near-neigbor Coulomb interactions are retained. The filling de- 
pendence of the d x 2_ y 2 transition temperature is studied in both the "hole-doped" and "electron- 
doped" regimes using parameters derived from constrained-occupancy density-functional theory for 
La2Cu04. The agreement with experiment on the overdoped hole side of the phase diagram is 
remarkably good, i.e., transitions emerge in the 40 K range with no free parameters. In addition the 
importance of the "orbital antiferromagnetic," or flux phase, charge density channel is emphasized 
for an understanding of the underdoped regime. 
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I. INTRODUCTION 



An experimental consensus has developed in recent years that the order parameter in the high-temperature cuprate 
superconductors has d x 2_ y 2 symmetry.EI Well before experiments indicated this esc etic symmetry a variety of theoretical 
approaches had suggested a tendency toward d x 2_ y 2 pairing in the HubbardoQ and t — J models.cl Within weak- 
coupling approaches, which treat the Coulomb interaction as a perturbation to one-electron band theory, exchange of 
antiferromagnetic spin fluctuationsO leads to pairing. 

While the correctness of the spin fluctuation scenario remains controversial, it is of interest to examine the pairing 
process within a more realistic setting than the one-band Hubbard model. It is well-established that magnetism in the. 
"undoped" cuprates can be understood within the context of a three-band modelQ (which projects to a t — J modelu 
in the strong-coupling limit). This CuC>2 model describes nearly filled Cu 3d x 2_ y 2, O 2p x , and O 2p y orbitak. which 
form a two-dimensional square Bravais lattice with a three-atom unit cell. The largest Coulomb integralsErEj in the 
CuC>2 model are the repulsion between holes on the same d orbital (Udd ~ 10 eV) or p orbital (U pp ~ 4 eV), and the 
repulsion between holes on neighboring d and p orbitals (f/ p d ~ 1 eV). 

A self-consistent and conserving calculation of one-particle properties ipnftc CuC>2 model based on exchange of 
magnetic and charge density fluctuations has been carried out previously£3li3 In the present paper we extend this 
fluctuation exchange, or "FLEX," calculation to an analysis of eigenvalues of the particle-particle and particle-hole 
vertex functions and the resulting transition temperatures. In particular, this analysis is carried out using one- and 
two-particle matrix elements deduced from constrained-occupancy density functional theoryB with no additional model 
projections or parameter fits. The results of this calculation with no adjustable parameters arc, if not compelling, at 
least suggestive. 

While the FLEX approach is inherently approximate, the observed trends in eigenvalues and transition temperatures 
for variations in hole density and Coulomb integrals can be expected to be carried over in more exact treatments. In 
addition this calculation provides a detailed example of the melding of many-body and band theory techniques now 
possible. 

The paper is organized as follows: The model and calculational notation are summarized in Section The particle- 
particle and particle-hole vertex functions within the FLEX approximation are derived in a computationally tractable 
form in Section [ill]. After a brief digression on sources of error, results for eigenvalues, transition temperatures, and 
eigenfunctions are presented in Section IV. The implications of the calculation are discussed, along with an overall 
summary, in Section 
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II. MODEL AND NOTATION 



In this section we define our notational conventions for the model to be studied. The three-orbital model for 
superconducting cuprate layers may be written in terms of creation operators for holes or electrons. As in Reference [l^, 
which we hereafter denote "EB," we adopt the hole representation; as an example, ct (R) creates a 3d x 2_ y 2 hole with 
spin ct in unit cell R. In addition we choose a staggered orbital phase which helps simplify the analysis of two-particle 
eigenstates. The unit cell and phase conventions are illustrated in Figure [j]. The Hamiltonian is conveniently broken 





Cu d x 2_y2 



CO ?* 



Op, 



FIG. 1. Unit cell and orbital phase conventions. The unit cell contains three orbitals: the Cu 3d x 2_ y 2, the O 2p x on x-axis 
bonds, and the O 2p y on y-axis bonds. The orbital phases are chosen in a checkerboard pattern. This assures that near-neighbor 
Cu-O and O-O hopping integrals have the same sign in all unit cells and greatly simplifies the two-body eigenstate analysis. 



up into one-particle and two-particle components 

H = H + V . 

With our conventions the one-particle Hamiltonian takes the form 

Ho-fiN = 
(£d - jti) 



t 



n d (R) + (e d -fi + e)^2 jn x (R) 

R 

P dE [ c L( R )<w( R ) + c dff (R + £)c XCT (R) 



R 



4(R) V (R) + 4(R + y) V (R) + H.C. 
*ppH [4-( R ) c m ( R ) + c yff (R + x)c X(T (R) 



a, R 



(1) 



n y (R) 



+ 4-( R -£) c M ( R ) + ct CT (R + £-z/)c XCT (R) + H.C. 
The number operators n d , n x and n y are defined in the usual way, e.g., 

n d (R) = £ c L( R )c dff (R) • 



(2) 



(3) 



The physical values of the short-range hopping matrix elements £ pd and t pp are both positive. Erli-J The d-hole creation 
energy e d may be set to zero without loss, and the p-d energy level difference e is positive. 

In the two-body-.Hamiltonian V we retain the three largest Coulomb integrals from constrained-occupancy density 
functional studies aa The on-site Cu repulsion C/ dd , the on-site O repulsion U pp , and the near-neighbor Cu-0 repulsion 
U p d- The last interaction complicates the analysis since it has both intra-cell and inter-cell components. The Coulomb 
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interactions may be written in a spin-diagonalized form which allows a decoupling of S — (density) and 5 = 1 
(magnetic) excitations. This procedure is treated at length in EB. 

The one-particle propagators for the Hamiltonian described above take the form 

G ab (R a , t q ; R 6 , n) = G Qb (AR Qb , Ar ab ) = -(T r c a (R a , r a )4(R 6 , t 6 )) , (4) 

where (a, R a , r a ) and (b, R b , 7b) are the orbital, unit-cell, and imaginary-time labels for particles in the final and 
initial states (see Figure g). Reference |l^ describes the general procedure for calculating such propagators and provides 



b 

R b 

FIG. 2. Diagrammatic representation of the one-particle propagator GatfARnt, Ar a b)- 

detailed results for the CuC>2 model described above. For the remainder of this paper we make use of one-particle 
properties obtained in this previous study. 



R 



III. DERIVATION OF VERTEX FUNCTIONS 



The calculation of eigenvalues of the particle-particle kernel in the Cu02 model is conceptually straightforward, 
but notationally involved. It is assumed that self-consistent one-particle propagators G have been obtained using the 
technique described in EB. Functional differentiation of the off-diagonal self-energy in the presence of an external 
pairing field yields the irreducible particle-particle vertex r pp . Using the notation developed in Reference [li] the 
singlet and triplet parts of the vertex are as follows (see Figure |J) : 

r PP (12; 34) = 14(12; 34) 

+ i$ d (24; 31) - f $ m (24; 31) + ±$ d (14; 32) - |$ m (14; 32) , (5) 
r t pp (12; 34) = 1^(12; 34) 

+ i$ d (24; 31) + i$ m (24; 31) - ±$ d (14; 32) - |$ m (14; 32) . (6) 

The numerical indices represent the space and time degrees of freedom of each particle, i.e., 




+ {v d ,D^V m ,M} 

FIG. 3. Irreducible singlet vertex function Ff p within the FLEX approximation. Outgoing states are represented on the 
right of the diagrams, incoming states on the left. (The coefficients 1/2 and —3/2 are omitted for clarity; see Equation ffy.) V s 
is the unrenormalized Coulomb matrix element in the singlet channel. The vertical ladders represent the exchange of density 
fluctuations. 
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1 = (mi, Ri, n ) 



(7) 



with mi the orbital, Ri the unit-cell displacement, and t± the imaginary time coordinate for particle 1. 
The matrix functions <I>d and <I> m represent particle-hole ladders in the density and magnetic channels: 

$ r (12; 34) = [l/ r G ph (l - VrG^y 1 ^] (12; 34) 



(8) 



for r = d and m. The matrices V T are the spin-diagonalized Coulomb interactions in each channel, and the matrix 
G ph is the uncorrelated particle-hole propagator: 



G ph (12; 34) =/3(T T c(l) C t(3)) (T T c(4) C t(2)) 
= /3G(13)G(42) , 



with P the iaverse temperature. 

As usual ]13 matrix multiplication is defined by 



(AB) (12; 34) = A(12; 56) B(56; 34) 



(9) 



(10) 



with an implied sum on repeated indices. The singlet and triplet kernels are obtained by multiplying the vertex 
functions by the uncorrelated particle-particle propagator 



G pp (12; 34) 



±/3G(13)G(24) 



(11) 



Note the presence of ^ in this definition of the propagator, which is consistent with our normalization of the vertex 
functions below. 

Expressions for the density and magnetic Coulomb matrix elements Vd and V m have been given previously in EB. 
Explicit expressions for V s and Vt follow from the diagrams in Figure § As in our previous work, it is convenient 




FIG. 4. Representation of the unrenormalized singlet and triplet Coulomb matrix elements V s and Vt. 

to adopt a notation which emphasizes the dependence of the matrix elements on only three unit-cell displacements 
between the two initial-state and two final-state particles (see Figure ^|). Thus, V s (ATl ac ; ab, AR„i,; cd, AR C( j) is the 
singlet Coulomb matrix element for a final-state particle pair in orbitals a and b with relative unit-cell displacement 

AR ab = R a R b , (12) 

and an initial-state particle pair in orbitals c and d with relative unit-cell displacement AR C£ 2. The displacement 
between the initial and final state particles is given by AR ac . There are only 11 two-particle states (ab, AR a fc ) which 
have non-zero singlet and triplet Coulomb matrix elements in the Cu02 model considered here. These states are 
listed in Table [j| for the unit-cell depicted in Figure |l|. (An identically labeled 11-state basis for non-zero density and 
magnetic Coulomb matrix elements was defined in EB.) 

As in EB, the initial/final state displacement AR QC is conveniently eliminated in favor of a center-of-mass momentum 
Q by writing 
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FIG. 5. Definition of unit-cell displacements in the representation of the singlet Coulomb matrix element. 
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TABLE I. Indexing scheme for the minimum-range particle-particle basis set in the CuCh model. The particle orbitals are 
a and b, with corresponding unit-cell displacement AR a t, = R„ — Rt,. Note that kernel eigenstates must satisfy the symmetry 
requirements of the Pauli Principle, but the basis states need not. 



V S (Q; ab, AR afc ; cd, AR cd ) = ^ g — tQ AR oc T/ s(ARac . a6j ARafc . cdj AR C(i ) . (13) 

AR« 

The indices in Table | can then be used to write V s and 14 compactly as Q-dependent 11 x 11 matrices. For example, 

V S 33 (Q) = 2f/ pp 

K 64 (Q) = -vf 4 (Q) = t/pd 

Vj 5 (Q) = -V t 75 (Q) - e^U pd . (14) 

Though the basic Coulomb interactions V s and Vt are short-ranged, the fluctuation-induced contributions to the 
particle-particle vertex functions T^ p and F^ p are not. Nevertheless, it is possible to calculate accurate pairing 
eigenvalues using vertex functions truncated in the relative displacement of the particle pair. For this reason it is 
convenient to arrive at a particle-particle vertex labeled using (i) total momentum-frequency Q = (Q, iQ); (ii) pair 
orbital indices (ab) (3x3 = 9 possible combinationsJbr the three-orbital model); (hi) unit-cell displacement AR a /, 
of the pair elements; and (iv) relative frequency ioj$3 (There is no additional benefit in introducing a relative time 
coordinate, since the fluctuations induce long-range couplings in imaginary time.) Previous notation for the time- 
independent Coulomb matrix elements may be generalized in a natural way. The desired singlet and triplet vertex 
functions (see Figure |^(a)) take the form 

rP p (Q; TOirn 2 , AR12, iuj; m 3 m 4 , AR 34 , iw') . (15) 



In order to calculate the crossed-channel particle- hole ladders $d and $ m , it is essential to use a different basis set 
obtained by a series of Fourier transforms. An initial Fourier transform on the relative displacement coordinates in 
Equation (||) yields 

r PP (Q; mim 2 , AR i2 , ilu; m 3 m,4, AR 34 , iuj') = 
Ks(Q; mxm 2 , ARi 2 ; m 3 m 4 , AR 34 ) 
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FIG. 6. Calculation of the irreducible particle-particle vertex functions Tf p , r = s and t. (a) Diagrammatic representation 
of the irreducible vertex in the computationally optimal basis set. Note that the total center-of-mass momentum-frequency 
Q = (Q, iQ.) is conserved, (b) Fourier-transformed singlet vertex function AF!? P ((2; m\m,2, k; 77137714, k'). See also Equa- 
tion (Q. (c) Representation of the first particle-hole ladder in (b) in the relative displacement basis. See also Equation (|is|). 



G 



eJ k-AR 12 Ar PP (Q; mim2) fe . m3m4) fe ') e 



/\ -ik'-AR 34 



(16) 



kk' 



with 



k = (k, iu>) 
k! = (k', io; ) . 

The fluctuation-induced contribution Ar PP takes the form 

Ar PP (Q; mim2, k; m^m^, k') — 

(k' — fc; m 2 m 4 , — k ; m 3 m!, k + Q) 

(k + k' + Q; TO1TO4, — fc'; 77137712, — fc) 



(17) 



2 d 2 m 



(18) 



The ladders are represented diagrammatically in Figure 6|(b). 

The first particle-hole ladder in Equation (|l8|) may be translated back to the relative displacement basis (Fig- 
ure H(c)): 

^d{k' — k; m 2 m 4 , — k ; 777 3 toi, fc + Q) = 



E 



ik'-AR^ 



$ d (fc' - fc; to 2 ™ 4 , AR 24 ; mam!, AR 31 ) e 4 ( k +Q)- AR 3i 



(19) 



Note that primes are included on the displacements here to emphasize that they are dummy summation variables, 
at this stage unrelated to AR42 and AR34 in Equation (jl^). Similar expressions hold for the other ladder sum 
terms contributing to Ar PP . Note that <£>d and $ m are independent of the relative frequency variables (due to the 
instantaneous character of Vd and V m ) and have been dropped from the notation without loss. 

The expressions in Equations (||) and (|l|) each involve double Fourier transforms and are impractical to calculate 
numerically. A much simpler form for Ar may be derived by changing momentum variables and interchanging the 
order of sums. For example, for the first particle-hole ladder contribution the appropriate change of variables is 

k' k -> Q' 

k' k' . (20) 

The sum on k' may then be carried out explicitly, yielding a delta function, which collapses the sum on AR 24 . 
After additional relabeling of summation variables, the complete result for T PP which results from this procedure is 

r PP (Q; TO1TO2, AR12, iu; m^m A , AR 34 , iu;') = 
V S (Q; mim 2 , AR12; m 3 m 4 , AR 34 ) 
+ e^ QARl3 Ar PP (ARi 3 , ifl; mim 2 , AR 12 , iu; m 3 ™4, AR 34 , iu') , 

AR 13 

(21) 

with 



Ar PP (ARi 3 , iQ; mim-2, AR12, iu; m^m^, AR34, iu 
1 



N 



E e iQ '' 



AR 2 



Q' 



J_ ST p iQ' AR13 

N ^ 



3*d ~ §*„ 



(Q', i(u' - u); m 2 m 4 , AR 24 ; m 3 mi, AR 3i ) + 
(Q', i(u + J + SI); mim 4 , AR i4 ; m 3 m 2 , AR 32 ) , 



Q' 



(22) 



where all relative displacements are expressed in terms of AR13, AR12, and AR3 4 : 
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AR.23 — — AR32 — AR13 — AR-12 

AR 24 = ARi 3 - ARi 2 + AR 34 
AR31 = -AR 13 

AR14 = AR13 + AR34 . (23) 

A similar expression for p may be obtained immediately using the correspondence in Equation (g) . 

The ladder summations and $ m may be calculated as matrix products in the space with compound indices 
(ab, AR ab ): 

$r = -VrXil + V^)- 1 ^ , (24) 

for r=d, m, where the uncorrected fluctuation propagator \ is defined by 

X(Q; ab, AR afc ; cd, AR cd ) = ~ £ e ^(AR ab -AR cd ) G ^ k + Q)Gdb (k) . (25) 

k 

For the Cu02 model the required matrix inverse is only If x ff. Note, however, that a separate inverse must be 
calculated for each value of the particle- hole ladder's center-of-mass momentum- frequency. 

The uncorrelated particle-particle propagator may also be expressed in the basis adopted above (Figure ^) : 



G PP (Q; mim 2 , ARi 2 , iw; m 3 m,4, AR34, iui' 
T 
N 

The particle-particle eigenvalue problem then takes the form 



-¥»»'Jx E e 4k - (ARl2 - AR34) G mim3 (k + Q)G m2rni {-k) . (26) 

k 



Q 

m, * m, 



i'"3 ' 1 i 

AR 34 t t AR 12 

I m, >. I 



-1(0 -ZOO 

FIG. 7. Diagrammatic representation of the uncorrelated particle-particle propagator 

G PP (Q; mim2, AR12, iu>; 77737714, AR34, iu'). Note the propagator is diagonal in the relative frequency, i.e., it vanishes for 



rr(Q)GPP(Q)0(Q) = \(Q)4>(Q) , (27) 

for r = s and t. Note that with the conventions adopted here a positive eigenvalue indicates attraction. (Although 
the kernel is non-hermitian, it is possible to show for Q — that the particle-particle eigenvalues are real-valued or 
occur in complex conjugate pairs.) The matrices T PP and G pp operate in a far larger compound- index space than that 
defined previously for the Coulomb interactions V T . The index now consists of the orbital-pair label (77117712), which 
takes on nine values in the CuC>2 problem; the subset of unit-cell displacements ARi 2 retained; and the set of values 
of the relative frequency to within a pre-defined cutoff interval. 

Note that since the kernel matrix is non-hermitian, its sets of left and right eigenvectors are not simply related 
(even though the left and right eigenvalue spectra are identical.) In the following section we emphasize the real-space 
and frequency dependence of the right eigenvectors, i.e., those determined by Equation (p7[). This is natural since the 
right eigenvector at T c evolves smoothly into the off-diagonal self-energy below T c . (The right eigenvalue equation 
may be re-derived by linearizing a self-consistent field problem in the off-diagonal self-energy) The corresponding left 
eigenvector has no such simple physical interpretation. 

A number of powerful approaches have been developed in recent years to compute a few selected eigenvalues of a gen- 
eral non-hermitian matrix in cases such as this for which a full diagonalization is impractical. All such approaches are 
derived from the much more standard algorithms available for the real-symmetric and complex-hermitian eigenvalue 
problems. We have made use of a so-called Lanczos-Arnoldi algorithm developed in the Department of Computational 
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and Applied Mathematics at Rice University.0 Using this algorithm we have studied kernels with row dimensions of 
order 10,000. 

To supplement our study of particle-particle eigenvalues we have also calculated a set of kernel eigenvalues for the 
particle- hole channels. These channels describe scattering of S = (charge density) and 5=1 (magnetic) excitations. 
Several points are important to note in this regard. First of all, the FLEX calculation (and any Baym-Kadanoff 
approachED) lacks self-consistency at the two-particle level. For this reason the density and magnetic propagators 
which enter the one-particle self-energy are not the same as those obtained by functional differentiation of the self- 
energy with respect to an external field. The difference may be described in terms of "vertex corrections" to the bare 
density and magnetic matrix elements Vd and V m . Within FLEX the simplest vertex corrections AFd and Ar m have 
a form closely related to the singlet and triplet interactions AT S and Ar t , i.e., they represent the exchange of single 
crossed-chaBnei density and magnetic fluctuations. More complicated vertex corrections take the Aslamazov-Larkin 
(AL) forrrdi3iI3 i.e., they describe the emission and re-absorption of pairs of fluctuations. For reasons described 
previouslyEj we omit the AL corrections to AFd and Ar m in the analysis which follows. 

It is also important to re-emphasize at this point that the one-particle FLEX calculations described here and in 
EB assume the exchange of elementary particle-hole fluctuations, but not elementary particle-particle fluctuations. 
For this reason particle-particle fluctuation propagators do not appear in crossed-channel contributions to AFd and 
Ar m below. In analogy with Equations (||) and (||) the spin-diagonalized particle-hole vertices (see Figure ||) may be 
written as follows: 



rP h (12; 34) = Vd(12; 34) - ±$ d (42; 31) - |$ m (42; 31) 
^(12; 34) = y m (12; 34) - ±$ d (42; 31) + ±$ m (42; 31) 

The functions <I>d an d 'I'm are as defined previously. 



(28) 
(29) 



+ 




+ 




FIG. 8. Irreducible density vertex function Fjj h within the FLEX approximation. Note the absence of AL and particle-particle 
exchange diagrams discussed in the text. (As before, the coefficients —1/2 and —3/2 are omitted for clarity; see Equation (pq).) 

In terms of the center-of-mass momentum-frequency Q the density vertex takes the form 

r^ h (Q; TOi77i 2 , AR12, iuj; m 3 m4, AR 34 , iui') = 
V d (Q; toito 2 , AR12; m 3 m 4 , AR 34 ) 
+ e- iQ ARl3 ArP h (AR 13 , mi m 2 , AR 12 , iui; m 3 m 4 , AR 34 , iui') , 



with 



AR13 



(30) 



ArP h (AR i 3 , ifl] m 1 m2, AR 12 , iui; m 3 m 4 , AR 34 , iui') 



N 



E 

Q' 



oi Q'-AR 4 



■i*d-§*„ 



(Q', i(ui' - u>); m 4 m 2 , AR 42 ; m 3 mi, AR 3i ) 



(31) 



where, as in Equation (p3[), all relative displacements are expressed in terms of the set {ARi 3 , AR12, AR 34 }. The 
analogous expression for r^ h follows by the correspondence in Equation 
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The particle-hole eigenvalue problem takes the form 

Tf(Q)G» h (Q)cf>(Q) = \(QMQ) , (32) 

where now 

G ph (Q; mim2, AR12, iuj; m 3 m4, AR34, iu) 1 ) = 

<W ^ e lk -( ARl2 - AR -) G mim3 (k + Q)G mim2 (k) . (33) 

k 

As in Equation (|27|), a positive eigenvalue indicates attraction. 

IV. RESULTS AND DISCUSSION 
A. Sources of Systematic Error 

In this section we discuss the nature of the errors which arise in calculation of instability eigenvalues and transition 
temperatures. Four main sources of error arise in the eigenvalue calculations. These are accumulation of frequency- 
space renormalization groupE2l error at low temperatures; error from the use of frequency cutoffs; error from truncation 
of the two-body vertex function in the relative real-space coordinate; and k-space discretization error (16 x 16 meshes 
are employed throughout). Detailed discussions of the renormalization group procedure for the one-particle self-energy 
are included in EB and Reference |2^. The errors associated with this approximation are generally negligibly small in 
comparison with the other sources. 

The frequency cutoff used in our calculations is fl c = 0.5£ p d for the ingoing and outgoing frequencies u> and u>' 
(see Figure ^) in the fluctuation-induced component of the singlet kernel. For the instantaneous part of the singlet 
kernel, whose decrease at high frequencies is controlled solely by the falloff of the uncorrelated propagator G pp , 
the corresponding cutoff is 50i p d- Errors associated with these cutoffs are extremely small. For example, for the 
standard parameter set (see Equation (|3~i|)) at 16% hole doping and temperature T — £ p d/512 (29 K), the d x 2_ y 2 
eigenvalue obtained using the cutoffs described above is Ad = 1.0458. If both cutoffs are raised to 50i p d, the eigenvalue 
becomes 1.0459, a change of 0.01%; this demonstrates the calculation's insensitivity to the cutoff associated with the 
fluctuation component. In contrast, if both cutoffs are dropped to 0.5i p d, the eigenvalue becomes 1.0439, a change of 
0.2%; this demonstrates insensitivity to the cutoff associated with the instantaneous component. It should be noted 
that at higher temperatures the cutoff on the fluctuation component must be raised to obtain comparable percentage 
accuracy. This is not costly, however, since the density of Matsubara frequencies decreases at the same time. 

Next we discuss the truncation procedure for dealing with the relative real-space coordinate in the two-body vertex. 
When the kernels are evaluated on a 16 x 16 k-space grid, the relative displacements AR12 and AR34 (see Figure |fj|) ma y 
take on 256 different values. Since the d x 2_ y 2 eigenfunctions fall off rapidly at large values of AR12 (see Section ^V D| ) , 
it is rather intuitive to introduce a truncated basis set for the relative displacements. In our calculations we limit the 
basis set to the twenty-one smallest lattice vectors; i.e., elements of the kernel are zeroed out for |AR| > a\/5. The 
corresponding gain in computation time is approximately (256/21) 2 ~ 150. 

Since the calculation of the full model with the untruncated real-space basis set is too time-consuming to be 
practical, we have used the simpler model with U pp = U p d = for an error analysis. The behavior of the two models 
is expected to be identical as far as this error check is concerned. In Figure || we plot the temperature dependence 
of the d x 2_ y 2 eigenvalue for the C/dd~ on ly model using the untruncated basis set and the 21-state basis set. The 
difference in the eigenvalues is very small for the two cases. For example, at T = i p d/1024 (15 K), Ad = 1.0683 with 
the untruncated basis and Ad = 1.0572 with the 21-state basis. The corresponding T c values are 20.8 K and 20.1 K, 
justifying the use of the truncated basis set. 

The biggest source of error in the calculation of the instability eigenvalues is the use of a 16 x 16 k-space discretization. 
For the models under study, the low-temperature eigenvalues from a 16 x 16 and a 32 x 32 discretization differ bjjJesa 
than 5%. This discretization error is very similar to that in previous studies of the one-band Hubbard modelJlSEII 
This means one should also expect roughly the same size error (i.e., 5%) in comparing the 16 x 16 results to the 
fine- mesh limit. 

In the figure below we plot the temperature dependence of the d x 2_ y 2 eigenvalue for the Z7dd-only model using 
16 x 16 and 32 x 32 discretizations. (Essentially identical behavior is expected for the the full CuC>2 model.) At 
T = < pd /1024 (15 K), A d = 1.1034 for the 32 x 32 study and 1.0572 for the 16 x 16 study. For both cases a 21-state 
real-space basis truncation has been employed. The corresponding T c values are 24 K and 20 K, corresponding to an 
underestimation of T c by 4 K using the 16 x 16 discretization. 
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FIG. 9. Systematic error analysis for eigenvalue and T c calculations in the C/dd-only model. All parameters are at their 
standard values (Equation (^)) except that U pp = U p d = 0. The filling is (n) — 1.16. (a) Comparison of temperature-dependent 
d x 2_ y 2 eigenvalues calculated using a full basis of relative displacement states (solid line) and the 21-state basis with | AR| < ayE 
(crosses). The k-space mesh is 16 x 16. (b) Comparison of eigenvalues calculated using a 32 x 32 discretization (solid line) and 
a 16 x 16 discretization (crosses). The 21-state truncated basis is employed. 
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As mentioned in Section [II we have employed a Lanczos-Arnoldi algorithm^ to calculate the first few maximum- 
real-part eigenvalues in each scattering channel. This algorithm is especially powerful for sparse matrices because it 
requires only repetitive multiplication of a vector by the matrix of interest. Since a large fraction of the elements in 
our scattering kernels are non-zero, but negligibly small with regard to calculation of the large eigenvalues, a sparse 
storage scheme is appropriate. For the scheme adopted throughout most of our calculations, eigenvalues are affected 
by less than a few parts in a thousand, and the gain in storage is of order 50. 



B. Eigenvalues for Particle-Particle Channels 

In the plots which follow we Make use of a "standard" CuC>2 parameter set derived for undoped La2CuC>4 by Hybert- 
sen, Schluter, and ChristcnsenEI using constrained-occupancy density functional theory. These standard parameters 
for the Hamiltonian in Equations (|l|) and (||) are as follows: 

i pd ~ 1.3 eV = 15, 100 K 

t pp ~ 0.65 eV = 0.5< pd 
e ~ 3.6 eV = 2.75i pd 
U dd ~ 10.5 eV = 8t pd 
£/ pp ~4eV = 3i pd 

U pd ~ 1.2 eV = t pd . (34) 

The temperature dependence of the maximal particle-particle eigenvalues for the standard parameter set at (n) = 
1.16 (16% hole doping) is illustrated in Figure [HJ. The maximal singlet eigenvalue corresponds to a d x 2_ y 2 state. 
This eigenvalue reaches unity, indicating a superconducting transition, at T/t pd = 0.0025, i.e., T — 37 K. At the 
transition temperature the next-leading singlet eigenvalue is of order 0.4 and corresponds to a state with so-called 
g-wave symmetry (i.e., nodes on the x and y axes, as well as the lines x — ±y; see Figure |TT| ) . A third eigenvalue, 
corresponding to an orthogonal d x 2_ y 2 state, lies just below the g-wave. 




FIG. 10. Temperature dependence of the maximal singlet and triplet eigenvalues for the standard parameter set at (n) = 1.16. 
The singlet eigenfunction has d x 2_ y 2 symmetry, and the triplet state odd-frequency s-wave symmetry. The d x 2_ y 2 eigenvalue 
reaches unity, signaling a superconducting transition, at T = 37 K. 

The maximal triplet channel eigenvalue in Figure [l0| remains small (~ 0.2) throughout the temperature range of 
interest. The triplet state in this case is antisymmetric in frequency and s-wave- like (i.e., symmetric) in space. (Note 
in this regard that our instability analysis includes all eigenvectors of the scattering kernels, including exotic singlet 
and triplet states with an antisymmetric frequency dependence.) 
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xy(x 2 -y~) 




FIG. 11. Schematic representation of the nodal structure of the g xy ( x 2_ y 2j singlet state. Note that the subscript, when 
viewed as a function, vanishes on the locus of nodes (just as in the case of the d x 2_ y 2 state). 

For comparison the behavior of the maximal particle-particle eigenvalues at (n) = 1.00 is illustrated in Figure |l^. 
The extreme singularity of the magnetic fluctuations in this case prevents study at temperatures lower than i p( j/64, 
i.e., T = 240 K. 
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FIG. 12. Temperature dependence of the maximal singlet and triplet eigenvalues for the standard parameter set at (n) = 1.00. 
The eigenfunction symmetries are as in Figure NOl 



C. Transition Temperatures for d x 2_ y 2 Superconductivity 



Eigenvalue plots of the type illustrated in Figure |l0| may be used to extract transition temperatures for the d x 2 _ y 2 
singlet. The critical behavior of this FLEX transition is classical, despite the fact that it is driven entirely by 
fluctuations. In terms of the FLEX eigenvalues, 
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A d (T) ~ l-a{T-T c ) 



(35) 



for T r*> T c , with a > 0. This contrasts with the exact critical behavior for a two-dimensional superconducting 
transition in the xy universality class: 



, xy (T) ~ 1 - B(T)e- A /V^^ , (36) 

with A a positive constant and B(T) an algebraic function.il It is nevertheless possible to interpret the d x 2_ y 2 
instability in FLEX as a "mean-field" transition with respect to critical order parameter fluctuations. With this 
caveat it is of interest to examine the dependence of this instability on doping and model parameters. [A more 
sophisticated treatment of the interference between the d x 2_ y 2 transition and the incipient instability in the magnetic 
channel is presumably necessary for a detailed understanding of the pseudogap regime observed in experiments, elbut 
that is not ousdntention here. In fact an additional charge density state, the so-called "orbital antiferromagnet'x 



"flux phase, "E3 is also apparently relevant in the pseudogap regime; see the discussion of this state in Section IV E.] 
In the plots which follow transition temperatures are given in units of K; they may be rescaled in units of t p d using 
the correspondence in Equation (|34"|). The experimentally observed transition temperaturesc3 for La2- x Sr x CuC>4 are 
plotted for comparison using the assumed correspondence 



- 1 



(37) 



As shown in Figure [13], a d x 2_ y 2 transition occurs for both hole doping ((n) greater than 1) and electron doping 
((n) less than 1). Since the CuC>2 model has only approximate particle-hole symmetry around the point (n) = 1, the 
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FIG. 13. Doping dependence of the d x 2_ y 2 transition for the standard parameter set. Results are shown for both (n) >>-4..00 
(hole doping) and (n) < 1.00 electron doping. For comparison the doping dependence of the experimental transitionE£l in 
La2- x Sr x Cu04 is plotted (dashed line) using the assumed correspondence x — > (rt) — 1. The increase of the FLEX T c on the 
electron-doped side is largely due to an increase in the spin fluctuation strength. 

transition temperatures are not symmetric. Within our FLEX calculation the pairing interaction becomes increasingly 
singular as (n) — > 1, and we have only been able to calculate superconducting instability temp eratu res for doping levels 
greater than 12%. (In any case a self-consistent parquet-like treatment of vertex functionsllirllj seems essential for 
values of (n) closer to unity.) The higher transitions for electron doping are consistent with the presence of enhanced 
magnetic fluctuations on this side of the phase diagram. The transition temperatures on the hole-doped side are 
strikingly similar to the experimental curve in the overdoped regime, (n) — 1 > 0.16. At smaller doping the FLEX 
curve continues to rise, while the experimental curve peaks and turns down in the underdoped region. As remarked 
previously, in this region the d x 2_ y 2 singlet channel is in strong competition with the Q ~ (jr, tt) antifcrromagnetic 
spin channel, as well as an exotic Q ~ (it, tt) charge density channel (see also Section IV E). It is tempting to speculate 
that the downturn in the experimental d x 2_ y 2 transition temperature results from this competition. 
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FIG. 14. Dependence of T c on the 0-0 hopping integra^4 pp . In Figures |1Q| all model parameters are set to their standard 
values except as noted. Further the experimental curvec3 for doped lanthanum cuprate (dashed line) is superimposed for 
comparison. 



In the next six figures we examine the sensitivity of the d x 2_ y 2 transition temperature to changes in the model 
parameters. Our discussion is limited to the hole-doped side of the phase diagram. First we alter a single parameter at 
a time, keeping other parameters fixed at their standard values, then we briefly consider the behavior of the drastically 
simplified CuC>2 model with U pp — U p d = 0. 

The effect oLremoving the 0-0 hopping integral t pp is shown in Figure |l4|. This change alters the shape of the 
Fermi surface,E2l improving the degree of nesting and enhancing the spin fluctuation spectrum. However, the transition 
temperature remains essentially unchanged, since the positive effect on the singlet vertex is largely compensated by 
a reduction in the uncorrelated propagator G pp . 

The effect of changing the Cu-0 orbital separation e = e p — £d is much more drastic, as expected. The value of 
e largely determines the strength of thg spin fluctuations. (This is because e is smaller than C/dd, i-e., the system is 
in the so-called charge-transfer regime]]) For small values of e, occupation of the O orbitals becomes comparable to 
occupation of the Cu orbitals (or even larger, when Coulomb interactions are taken into account). As an example, 
for e = and (n) = 1.16, only 33% of the holes reside on the Cu orbitals. Since U pp is considerably less than Udd, 
increased O occupancy reduces the strength of the spin fluctuation propagator and weakens the pairing tendency. 
This fact is illustrated clearly in Figure [L5|. The d x 2_ y 2 transition temperature drops sharply when e is reduced from 
3.6 eV to 2.0 eV. The transition disappears completely when e is set to zero (i.e., the bare Cu and O orbitals become 
degenerate) ; this is due not only to the reduction of the effective Coulomb parameter, but also to the loss of nesting 
in the e = Fermi surface. 

The dependence of T c on the Coulomb parameters Udd, U pp , and U p d is relatively complex, since these parameters 
contribute to both the one- and two-body effective interactions. For the simpler one-band Hubbard modelEil (and 
for the Cu02 model with Udd only — see the discussion below), an increase of the Coulomb integral leads to a peak 
T c , then a gradual decrease at larger values. The origin of this behavior is a competition between the pairing vertex 
(which is enhanced by a large Coulomb interaction) and the uncorrelated propagator G pp (which is suppressed). In 
these simpler models the on-site Coulomb interaction does not directly suppress pairing, since the d x 2_ y 2 state has 
no on-site pairing component. While this remains true for Udd in the full Cu02 model, it is not necessarily true for 
U pp and U p d'- the d x 2_ y 2 pair wave function generally has on-site 0-0 and near-neighbor Cu-0 components, which 
are suppressed by the Coulomb integrals U pp and U p d- The importance of this direct effect depends on the a dmix ture 
of the relevant components in the d x 2_ y 2 pair state (see the discussion of the pair wave function in Section IV Y\ ). 

For the 0-0 Coulomb integral U pp , this direct suppression of pairing apparently dominates, i.e., an increase in U pp 
leads to more repulsion in the d x 2_ y 2 pair state and a reduced transition temperature (Figure [l6| ). As discussed in 
Section [VD below, the d x 2_ y 2 pair does have a non-zero on-site 0-0 component, consistent with the observed trend 
in T r . 
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FIG. 15. Dependence of T c on the unrenormalized Cu-O level separation e = £ p — ea- The behavior of the model for e — 
was also examined, but no transition occurs in this case. 
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FIG. 16. Dependence of T c on the intra-orbital O-O Coulomb integral U p 
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For the Cu-0 Coulomb integral U p d, the trend in T c with increasing U p d (Figure |T^) resembles the trend with 
increasing U in the one-band Hubbard model: An increase, peak, and subsequent decrease. This behavior reflects 




FIG. 17. Dependence of T c on the near-neighbor Cu-0 Coulomb integral U p d- 

a compromise between the enhancement of the spin fluctuations and pairing interaction with increased C/ pc j; the 
accompanying suppression of G pp ; and the direct suppression of T c noted above. The enhancement of the spin 
fluctuations with increasing £7 p d arises from improved Fermi surface nesting and from increased d-orbital occupancy 
due to a rise in the Hartree-Fock level separation e HF (see the discussion of the Hartree-Fock Hamiltonian in EB). 
[It is interesting to note that the direct suppression of T c is moderated by the sppe phenomenon which reduces the 
effective Coulomb repulsion in conventional electron-phonon superconductivity:EZI The near-neighbor component of 
the pair wave function changes sign at high frequencies, effectively reducing the repulsion in the low-frequency region, 
i.e., inducing a Coulomb pseudopotential.] 

Finally Figure |l8| shows the transition temperatures for two different values of Udd, 10.5 eV (standard parameter 
set — charge-transfer regime) and 2.5 eV (Hubbard regime-). The substantial reduction in Udd in this case has only 
a minimal effect on the strength of the spin fluctuations .EHl This is because a decrease in Udd results in an increased 
Hartree-Fock level separation £ H f and an increased d-orbital occupancy (cf. the discussion of an increase in U p d 
above). The increased d-orbital occupancy offsets the direct effect of a smaller Udd in the spin fluctuation propagator. 
The increase in T c for this admittedly unrealistic parameter set then results from an increase in G pp (i.e., a density 
of states effect). ._. 

The CuC-2 model with U pp — U p d = (the "C/dd-only model" ) has been studied previously^ at temperatures above 
the d x 2_ y 2 transition. This model is conceptually problematic: the omission of the interactions associated with the 
p-orbitals substantially alters the Hartree-Fock Fermi surface, largely negating any improvement in the band structure 
expected from the addition of the extra bands. (Both G pp and the d-orbital spin fluctuation strength are significantly 
affected by the omission.) Furthermore, while the d x 2_ y 2 wave function is dominated by d-orbital components, the 
omission of U pp and U p d completely eliminates those components associated with the p-orbitals. The principal virtue 
of the f/dd-only model is its calculational simplicity. Since the Coulomb interaction Udd is zero-range, the computations 
involved are essentially the same as those in the one-band Hubbard model. For example, the particle-hole ladders $ 



in Section III become scalar, rather than matrix, inverses. 

For completeness, the variation of T c with Udd in the C/dd-only model is shown in Figure [l9| As expected, the 
qualitative dependence of T c is the same as that in the one-band Hubbard modcl:Eil The peak value of T c for increasing 
Udd is determined by a competition between enhancement of the pairing interaction and suppression of the uncorrelated 
propagator G pp . 
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FIG. 18. Dependence of T c on the intra-orbital Cu-Cu Coulomb integral Um- 



D. d x 2_ y 2 Singlet Wave Function 



The graphical representation of the particle-particle pair eigenfunction 4>(Q; ab, AR a f,, iui) is hindered by its large 
number of degrees of freedom. It is essential to make use of symmetries to emphasize the eigenfunction's key features. 
The first basic symmetry follows from the Pauli Principle. Written in terms of anticommuting c-numbers, the pair 
state corresponding to eigenfunction </> r is 



e^ R+AR -> </> r (Q; ab, AR ah , iu) 

R; ab, AR ab , ui 



(38) 



where r=s or t, with 



and 



Since 



it follows that 



Xt = 



_L f 1 

" >/2 W 



i o\ j_ / o 1 

Oj' ^2 V 1 



or 




1 



Y, e jQ -( R+AR -' MQ; ab, AR afc , tw) 

R; ab, AR„i,u 

x Z Xr' a c fc< x(R, -iw)c aCT /(R + AR o6 , + = 

crcr' 

Z e^-( R+AR -) r (Q; a6, AR ah , iw) 

R; ab, AR ab , w 

x ^ x a / c a( x(R + AR ab , i(w + 0))c bCT / (R, -iw) . 



(39) 

(40) 
(41) 



(42) 
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FIG. 19. Dependence of T c on the intra-orbital Cu-Cu Coulomb integral Udd in the "t/dd-only" model. The Coulomb 
integrals U pp and U p d are in this case set to zero, with other parameters remaining at their standard values, (a) Variation of 
T c with doping for several values of Udd- The behavior for Udd/t p d = 1 was also examined, but no transition occurs in this 
case, (b) Variation of T c with Udd at fixed filling (n) = 1.12. 
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It is convenient to relabel the dummy sums on the left by first interchanging orbital indices a and b, then letting 

AR ba = -AR a j, 

R = R' + AR a ^ 

-iuj = i{uj' + fi) . (43) 
Dropping the primes on the dummy variables R' and u>' gives 

e ' Q ' R MQ\ ba > -AR ab , + n)) 

R.; ab, AR a (, , cj 

x 53 X^ V c aCT (R + AR ab , i(uj + n))c ba , (R, -iu) = 
- e jQ - (R+AR » b) r (g; ab, AR ab , iw) 

R; ah, AR a i,u 



J] xr'c aCT (R + AR Qb , £(w + n))c 6<r /(R, -iw) . (44) 



x 

aa 



The Pauli symmetry relations follow by identifying coefficients: 

MQ; ba, -AR ai , + ft))x r CT ' a = -e iQ - AR «^ r (Q; aft, AR ab , iw)^ , (45) 

i.e., 

UQ; ab, AR Qb , iw) = e- iQ,AH »»&(Q; 6a, -AR ab , -i(w + fi)) (46) 
for the singlet channel, and 

&(Q; aft, AR ab , %u) = -e-^ AR ^cj) t {Q; ba, -AR ab , + fi)) (47) 
for the triplet channel. These symmetry relations assume a particularly simple form for the case of interest here, 

= 0. 

To aid in the graphical display of the pair wave function it is useful to introduce a basis of one-particle states with 
simple transformation properties under point group operations. While the Cu 3d x 2_ y 2 orbital transforms into itself 
under all symmetry operations, the O 2p x and 2p y orbitals are generally mixed. It is, however, possible to form 
linear combinations of the p x and p y orbitals which transform in simple ways. To derive the transformed orbitals we 
rewrite the pair wave function 4> in Equation (|38|), holding the coordinates of particle a fixed. For brevity the spin 
and frequency dependence of <j> is temporarily suppressed. The components of the wave function for b — x and y take 
the form 

£ e lQ - (R+AR - ) c Q (R + AR ab ) 

R.; a, AR fl (, . uj 

x U r (Q; ax, AR ab )c K (R) + r (<2; ay, AR Qb )c„(R) . (48) 
In the first term the sums on R and AR ab may be shifted by 

AR ab -» AR Qb + x . (49) 

A similar shift may be performed in the second term with x — > y . This results in an equivalent symmetrized version 
of the wave function, 



i J- e^( R+AR -)c a (R + AR ab ) 

R.; a, AR a b. UJ 

Uh(Q; ax, AR ab )ca;(R) + r (Q; ax, AR ab + x)c x (R - x) 
r (g ; ay, AR afc )c y (R) + cj> r (Q\ ay, AR ab + y)c y (R - y)] (50) 



x 
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The four c- numbers c x (R) , c x (R — x), c y (R) , c y (R — y) may now be re-expressed in terms of the linear combinations 



[c D (R)- 




r 1/2 


1/2 


1/2 


1/2 " 




- c x (TL)^ 


c S (R) 




1/2 


1/2 


-1/2 


-1/2 




c x (H- x) 






1/V2 


-1/V2 










Cy(R) 


.Cy(R). 










-1/V2 


1/V2 




_c y (R-y) 



(51) 



The c-number 5d(R) represents an extended oxygen orbital with d x 2_ y 2 rotational symmetry (Figure pp|(a)), just 
like the central Cu 3d x 2_ y 2 orbital. (The uniform phases in the linear combination result from the initial definition of 
the 2p x and 2p y orbitals. ) Likewise c$ (R) represents an extended s-wave oxygen orbital (Figure pO|(b) ) , and cx (R) an d 
cy(R) represent extended p-wave orbitals (Figures ^(c) and (d)). The wave function components in Equation j50| ) 
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FIG. 20. Extended p-orbital basis set with well-defined rotational symmetry. The central Cu site is denoted by a shaded 
circle, (a) Extended d x 2_ y 2 orbital (state D). (b) Extended s orbital (state S). (c) Extended p x orbital (state X). (d) Extended 
p y orbital (state Y). 

may now be rewritten in terms of the new c-numbers. The complete pair wave function in Equation (|38|) then becomes 



£ e' Q - (R+AR " l) MQ; aB, AR Qb , iw) 

R; aB, AR,i,u 

XT*' Caa (R + AR afc , i((J + 0)) C Ba , (R, -ioj) , 



(52) 



where, as before, the sum on a runs over {d, x, y}, but now the sum on B runs over {d, D, 5*, X, Y}. The new wave 
function components are 

4> T (Q] aD, AR a6 , iuj) = 

4> T (Q; ax, AR afc , ioS) + ^ r (<9; ax, AR ab + x, iu) 

(j>r{Q; ay, AR a6 , iu) + r (Q; ay, AR a6 + y, iu>) 
(Q; aS, AR ab , iuj) = 

r (<3; ax, AR a fc, iu) + 4> T (Q; ax, AR ab + x, iuj) 
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4>r{Q] ay, AR Qb , iui) - 4> r {Q] ay, AR Qb + y, iui 
r (Q; aX, AR Qb , iui) = 

T (Q; ax, AR Qb , iui) - (j> T (Q; ax, AR ab + x, iui) 
(Q; aY, AR ab , iui) = 

h{Q\ ay, AR ab , iui) + (f> T (Q; ay, AR ab + y, iui) 



1 

2V2 



1 

2\/2 



(53) 



Note that this expression for the pair wave function is equivalent to that in Equation (|3q) . The new basis for the 
S-particles is simply overcomplete; the c-numbers cb ct '(R, iui) and cg CT '(R', iw) for near-neighbor R and R' are no 
longer independent. 

The new basis is well-adapted for representing pair wave functions for Q = in a simple way. It is convenient 
to keep the B-particle orbital and unit-cell position fixed while varying a and AR ab . For example, separate plots 
describe the system for the i?-particle in the Cu 3d x 2_ y 2 orbital {4> a d), in the extended oxygen d x 2_ y 2 orbital (4> a D), 
and in the extended oxygen s orbital (<f> a s)- If the pair wave function is to have overall d x 2_ y 2 symmetry, <\> a d and (j> a D 
must have explicit d x 2_ y 2 symmetry in the variable AR Qb , while <p a s must have explicit s symmetry. (States with d 
symmetry may also be constructed for 6-particles in the X and Y orbitals. These states are more complicated, since 
both (j) a x and 4> a Y must be non-zero.) 

In Figure [2l] we use histograms to show the spatial variation of the minimum-frequency (lo — irT) components of 
the d x 2_ y 2 singlet wave function (i.e., the right eigenvector of the particle-particle kernel) for T ~ T c . Recall that 
this function evolves smoothly into the off-diagonal pair field below T c . Each histogram shows a 4 x 4 patch of unit 
cells, centered on a Cu site in the cell at AR ab = 0. The orientation of the x and y axes is indicated in Figure pT](a) . 
The height of the block at each point in the lattice is just the value of 4> at that combination of a-particle orbital 
and displacement indices (a, AR ab ) for fixed i?-particle indices. It is clear that the wave function is dominated by 
the d-orbital components; the p-orbital components, however, play an important role in determining eigenvalues and 
transition temperatures, and their neglect is not justified. 

Finally the relative frequency dependence of several short-range components of the-, pair wave function 
4> s (0; ah, AR ab , iui) is displayed in Figure Note that, as in the one-band Hubbard model,Eil the wave function is 
strongly frequency dependent, falling rapidly to zero on a scale of ui ~ 0.5t p d- 



E. Eigenvalues and Wave Functions for Particle-Hole Channels 

A complete FLEX analysis of the particle-hole (i.e., magnetic and charge density) channels with the same degree of 
rigor applied in the particle-particle analysis is not attempted here. The reasons are as follows: (i) The version of FLEX 
considered in the present work and in EB is based on particle-hole exchange. Consequently the particle-particle vertex 
functions analyzed in the preceding section contain only single-fluctuation-exchange ladders, yet are fully conserving. 
On the other hand, a fully conserving calculation of the particle-hole vertex functions within this approximation 
scheme requires the inclusion of not only single-exchange ladders, but also a class of double-exchange Aslamazov- 
Larkin diagrams l 18 H 19 l Such diagrams constitute the beginning of a parquet-like renormalization of the single-exchange 
ladders .Ej Since this renormalization is incomplete (and does not improve the consistency of the particle-hole vertices 
which appear at different points in the calculation), the treatment of the Aslamazov-Larkin diagrams is problematic, 
(ii) In order to treat the particle-hole vertices on the same footing as the particle-particle, the FLEX approximation 
should include particle-particle exchange diagrams from the outset (see, e.g., References p^ , |lq , |l9| ) . Such a treatment, 
while in principle quite straightforward, exceeds the scope of the present work. 

The limitation imposed by these points makes a satisfactory analysis of the nearly singular magnetic channel 
impossible in the present work. This is because both the double- exchange Aslamazov-Larkin diagrams and the single- 
exchange diagrams from the crossed particle-particle channel are repulsive in the magnetic channel. The omission of 
these contributions in a naive calculation leads to a drastic overestimate of the magnetic eigenvalue (i.e., valuesdpger 
than unity). A similar situation has been noted in previous FLEX studies of the one-band Hubbard modeljuiia in 
that case the magnetic eigenvalue drops well below unity when the omitted contributions are reinstated. 

Note that these limitations in the treatment of the particle-hole channel do not compromise the conserving nature of 
the FLEX calculation for the particle-particle channel. (This is not to say that satisfying conservation laws guarantees 
accuracy: the overall lack of self-consistency in two-particle vertices which appear at different points in the FLEX 
calculation is a broader global concern Jiil which can be remedied only by a more sophisticated parquet-like analysis. 
See Section |V| for further comments on this point.) 
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FIG. 21. Relative displacement dependence of the minimum-frequency (w = ttT) components of the d x 2_ y 2 pair wave 
function at T = t pt j/512 = 29 K for the standard parameter set at (n) = 1.16. The histograms are centered on a Cu site 
in the cell with AR = 0. Note that the wave function is real-valued, (a) Orientation of the x and y axes in the histogram 
plots, (b) Component <f> s (0; ad, AR, inT). (AR labels the unit-cell displacement to orbital a from a fixed Cu 3d x 2_ y 2 orbital.) 
Note that the histogram exhibits explicit d x 2_ y 2 symmetry in this case, (c) Component 4> s (0; aD, AR, inT). In this case the 
fixed orbital is the extended O 2p linear combination with d x 2_ y 2 symmetry. As before the histogram exhibits explicit d x 2_ y 2 
symmetry. Note the difference in scale between this plot and that in (b). (d) Component 4> B (0; aS, AR, iirT). In this case the 
fixed orbital is the extended O 2p linear combination with s-wave symmetry. In this case the histogram exhibits explicit s-wave 
symmetry (see the discussion in the text). 
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Figure ||l] (continued) 
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FIG. 22. Relative frequency dependence of short-range components of the d x 2_ y 2 pair wave function for the standard 
parameter set at (n) = 1.16. The dominant component, associated with pairing on near-neighbor Cu 3d x 2_ y 2 sites, drops to 
zero over a range determined by the spin fluctuations. 



Despite these caveats, we feel it important to emphasize in this, section a feature of the physics in the charge 
density channel which has received little attention in recent years.Ea The Q = d x 2_ y 2 state in the singlet channel 
has an analog at Q ~ (tt, tt) in the density channel. The presenca-qf such an analog or partner state is familiar 
in a simpler context: in the one-band negative-?/ Hubbard model,[3Ej the physics near half filling is dominated by 
a Q = s-wave singlet state and a Q = (tt, tt) charge density state. At half-hlling these states become exactly 
degenerate, constituting the components of a Heisenberg-like order parameter. The instability in both the singlet 
and density channels is driven by the attractive unrenormalized vertex —U. In the positive-L/ Hubbard model and 
the CuC>2 model, an analogous pair of potential instabilities is driven by the exchange of spin fluctuations. In this 
case the singlet state of interest has d x 2 _ y 2 symmetry. The partner state, which becomes exactly degenerate with the 
d x 2_ y 2 singlet at half-filling in the positive-?/ Hubbard model, is a Q = (tt, tt) charge density state which shares the 
discrete (La . y 2 rotational symmetry. This state has been previously considered in both weak- and strong-coupling 
studics.EjcJ We follow SchulzEJ in denoting this state an "orbital antiferromagnet" ; the name is natural since the 
state describes microscopic currents which flow around elementary plaquettes in the square lattice, with the direction 
of current flow staggered between adjacent plaquettes (see Figure]^). In strong coupling the corresponding stateEj'ES 
has been denoted a "flux phase." 




FIG. 23. Representation of the circulating charge density currents in an ordered orbital antiferromagnetic state. The O sites 
are omitted for clarity. 



Away from half-filling in the Hubbard model and at arbitrary fillings in the CuC>2 model, the exact symmetry 
between the d x 2_ y 2 singlet and the orbital antiferromagnet is broken. Furthermore, with the loss of perfect nesting in 
the band structure, the wave vector Q for the optimal charge density state becomes incommensurate. For example, 
for the standard parameter set at (n) = 1.16 the optimal Q vector is approximately (1, 0.875)7r. 

We have studied the temperature variation of the orbital antiferromagnetic eigenvalue within an inherently limited 
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approximation: namely, we have included in the charge density vertex only single-exchange diagrams describing 
magnetic and density fluctuations. As noted above in the comments on the magnetic vertex, this approximation is 
not entirely satisfactory, since both double-exchange Aslamazov-Larkin diagrams and single-exchange particle-particle 
ladders are omitted; however, this approximation does preserve the crucial feature which determines both the d x 2_ y 2 
singlet and orbital antifcrromagnetic eigenvalues, i.e., the exchange of nearly antifcrromagnetic spin fluctuations. 
The temperature variation of the d x 2_ v 2 singlet and optimal orbital antifcrromagnetic eigenvalues for the standard 



parameter set is shown in Figure 24 While A af is smaller than Ad, the eigenvalues remain very close down to the 
d x 2_ y 2 transition. From this analysis it becomes clear that a fully satisfactory treatment of the model, particularly 





FIG. 24. Temperature dependence of the d x 2_ y 2 singlet and orbital antiferromagnetic (OAF) eigenvalues for the standard 
parameter set. (a) Results for filling (n) = 1.16. The optimal wave vector for the OAF state is in this case Q = (1, 0.875)-7r. 
(b) Results for filling (n) = 1.00. The optimal wave vector is in this case Q = (n, it). 

in the underdoped regime, must describe the competition between the magnetic, d x 2_ y 2 singlet, and orbital anti- 
ferromagnetic channels. While we have no evidence that the orbital antiferromagnet is ever actually the dominant 
instability, it is tempting to speculate on its relevance, at least in conjunction with the d x 2_ y 2 singlet, for a description 
of the anisotropic pseudogap observed in many experimental studiesEa 
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FIG. 25. Relative displacement dependence of the dominant minimum-frequency {ui — nT) component of the Q = (tt, n) 
orbital antiferromagnetic wave function, <^oaf(Q; ad, AR, inT) for the standard parameter set at (n) = 1.16. The temperature 
is T = tpd/512. As in Figure ^l], the histograms are centered on a Cu site in the cell with AR = 0. In this case the wave 
function is complex-valued, even though the eigenvalue is real. Note that a and AR vary with b — d held fixed. Note also the 
fact that d x 2_ y 2 rotational symmetry is manifest for a — d and arbitrary AR, but that the rotational symmetry is hidden for 
a = x and y, as discussed at length in the text, (a) Real part of the wave function, (b) Imaginary part of the wave function. 



In Figure ^ we show the spatial variation of the real and imaginary parts of the Q = (tt, tt) orbital antiferromagnetic 
wave function 4> OA f{Q', ad, AR a j, icu = ittT) for the standard parameter set at (n) — 1.16 and T = i p( j/512. The 
d x 2_ y 2 rotational symmetry of the wave function is manifest for the components with a = d, but is hidden for the 
components with a = x and y. This is true for the following reason: The total wave function takes the form 

J- e^-( R+AR -)0 OAF (Q; ab, AR ab , %u) 

R; ab, AR„i,u 

X H xT'caa{R + AR a6 , i{cj + n))c ba ,(R, ILO) , (54) 

era' 

with 

1 (\ 



Q = (tt, tt) 

a = . (55) 

It is convenient to adopt the shorthand a for the rotated image of orbital a and R~ for the rotated image of R a , 
the unit-cell location of orbital a. When the wave function coordinates (unit-cell and orbital labels) are rotated, it 
is guaranteed that the compound d-orbital label (a, R a ) = (d, R) maps to (a, R~) = (d, R), where R is the rotated 
image of R. However, under successive rotations, the p x -orbital label (x, R) maps to (y, R) (rotation through tt/2); 
(x, R — x) (rotation through tt); and (y, R — y ) (rotation through 3n/2). A similar set of transformations holds for 
the p v -orbital label. Discrete d x 2_ y 2 symmetry requires that 
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Figure (continued) 

e Q R ^ 4>oaf(Q; a&, R~ - R~, uo) = 

e 2lB e jQ R " O af(Q; ab, R a - R & , icu) , (56) 

where 

R b = R 

R a = R + AR ao 

6 = 0, tt/2, tt, 3tt/2 . (57) 
When the phase factors e » and e lC ^' Ra are equal, the d x 2_ y 2 symmetry is manifest 

0oaf(Q; a&, R~- Rp iw) = e 2ie O af(Q; a&, Ra - Rb, iw) ■ (58) 

However, when 

e jQ ' 5 7 = -e^- 11 - , (59) 

the d x 2_ y 2 symmetry is hidden, i.e., 

O af(<3; ab, R~- Rp iw) = -e 2ie ^ af(Q; a&, Ra - Rb, ia>) . (60) 

This accounts for the seemingly anomalous transformation properties of the tftxd and <p y d components of the wave 
function in Figure 

The frequency dependence of several short-range components of the Q = (tt, tt) orbital antiferromagnetic wave 
function is illustrated in Figure ^6[ Note that the wave function is in this case intrinsically complex, although the 
eigenvalue is real. While the Pauli Principle does not dictate the transformation properties under lo — > —lu in this 
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FIG. 26. Relative frequency dependence of short-range components of the Q = (n, n) orbital antiferromagnetic wave function 
for the standard parameter set at (n) = 1.16. As in Figure G2J, the component associated with near-neighbor Cu 3d x 2_ y 2 sites 
drops to zero over a range determined by the spin fluctuations, (a) Real part of the wave function, (b) Imaginary part of the 
wave function. 
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case, the overall phase of the wave function for Q = (tt, tt) and ft = may still be chosen such that 

0oaf(Q; ab, AR, — iuj)* = cf> OAF (Q; ab, AR, iuj) . (61) 

This follows from a basic symmetry of the eigenvalue problem for Kf h = rP h G ph in Equation (|3^), viz., 

Kf h (Q; mim 2 , AR i2 , iuj; m^m^, AR 34 , iu')* = 

Kf (-Q; mim 2 , ARi 2 , — iuj] m^m A , AR 34 , — iu') . (62) 

This implies that if <j){Q; ab, AR, iu>) is an eigenfunction with eigenvalue X(Q) for total momentum-frequency Q, then 
4>{Q', ab, AR, — ilu)* is an eigenfunction with eigenvalue X(Q)* for total momentum-frequency —Q. For the case of 
interest here Q — —Q, the eigenfunction 4>(Q; ab, AR, iu>) is non-degenerate, and the eigenvalue is real. Thus, 

4>(Q] ab, AR, -iuo)* = a(f>{Q; ab, AR, iuj) , (63) 

with a a complex constant. The constant may always be set equal to unity by a simple rescaling of (f>, leading to 
the symmetry relation in Equation (|6l|). Symmetries for general Q near {tt, tt) may be examined by an extension 
of this argument. Note finally the overall similarity of the d x 2_ y 2 singlet wave function plotted in Figures |2l] [2^ 
and the orbital antiferromagnetic wave function plotted in Figures |25|-|26], with respect to both spatial and frequency 
dependence. 



V. SUMMARY 



Our results demonstrate that the fluctuation exchange approximation provides reasonable results for both the 
magnitude and doping dependence of the d x 2_ y 2 transition temperature in the overdoped regime, (n) — 1 > 0.16. While 
the level of quantitative agreement between FLEX and experiments is almost certainly fortuitous, it is important 
to emphasize several points in this regard: (i) For a wide range of model parameters clustered around the standard 
LDA set,EI FLEX transition temperatures are predicted in the range of 10 to 100 K. It is by no means obvious that 
this should be so, i.e., one might have imagined obtaining a range of transitions stretching over several orders of 
magnitude, (ii) The continued rise .of the FLEX d x 2_ y 2 eigenvalue for values of (n) approaching unity is consistent 
with previous Monte Carlo studies ,£3 which have demonstrated enhanced d x 2_ y 2 correlations even in, regions where 
long-range magnetic order is being established. It appears clear that a more sophisticated approachtj is essential to 
resolve the competition between the incipient instabilities in the antiferromagnetic spin, d x 2_ y 2 singlet, and orbital 
antiferromagnetic channels in the underdoped regime (n) — > 1.0. 

As emphasized in Section IV, the presence of large eigenvalues in the orbital antiferromagnetic channel is an 
unambiguous result of our analysis, despite the technical limitations of our approach for the particle-hole channels. 
It is important to note that the orbital antiferromagnetic channel becomes degenerate with the d x 2_ y 2 singlet to 
form a Heisenberg-like order parameter in models with exact particle-hole symmetry (such as the half-filled one-band 
Hubbard model). While the breaking of particle- hole symmetry in the standard Cu02 model apparently favors the 
d x 2_ y 2 singlet at half- filling, it seems clear that the orbital antiferromagnet must be retained in any analysis which 
aims at a quantitative description of the region near (n) = 1 . 

Finally, as a more general comment, the present study demonstrates the feasibility of extending the FLEX instability 
analysis to models with an increasing degree of realism. While the principal shortcoming of FLEX, viz., the lack of self- 
consistency at the two-body level, remains a separate concern, it is also clear that progress toward a truly predictive 
many-body theory demands the ability to incorporate realistic details of lattice and interaction structure. A natural 
next step in this direction is the analysis of a one-band model with longer-range interactions. The general formalism 
developed in the present work and in EB (in particular, the use of a real-space basis set for relative coordinates) 
provides a calculationally feasible framework for such a study. 
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